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The computational study of chemical reactions in complex, wet environments is critical for ap¬ 
plications in many fields. It is often essential to study chemical reactions in the presence of applied 
electrochemical potentials, taking into account the non-trivial electrostatic screening coming from 
the solvent and the electrolytes. As a consequence the electrostatic potential has to be found by 
solving the generalized Poisson and the Poisson-Boltzmann equation for neutral and ionic solutions, 
respectively. In the present work solvers for both problems have been developed. A precondi¬ 
tioned conjugate gradient method has been implemented to the generalized Poisson equation and 
the linear regime of the Poisson-Boltzmann, allowing to solve iteratively the minimization problem 
with some ten iterations of a ordinary Poisson equation solver. In addition, a self-consistent pro¬ 
cedure enables us to solve the non-linear Poisson-Boltzmann problem. Both solvers exhibit very 
high accuracy and parallel efficiency, and allow for the treatment of different boundary conditions, 
as for example surface systems. The solver has been integrated into the BigDFT and Quantum- 
ESPRESSO electronic-structure packages and will be released as an independent program, suitable 
for integration in other codes. 


I. INTRODUCTION 

Many important chemical processes take place in soln- 
tion both in the context of basic and industrial research. 
The computational study of such chemical reactions in 
wet environments is therefore of cross-disciplinary inter¬ 
est to physics, chemistry, materials science, chemical en¬ 
gineering, and biologyh Gomputational studies can com¬ 
plement these investigations by giving insight into new 
processes and materials as well as reducing development 
times and production costs. Solar-energy harvesting in a 
dye-sensitized cell or electro-catalytic water splitting are 
two simple examples of relevance for applications in the 
energy and environment context. 

Molecular properties in the presence of a solution are 
often very different compared to pure in vacuum^ con¬ 
ditions making vacuum-like ab initio calculations an in¬ 
appropriate approach for such problems. An inclusion of 
the solute-solvent interaction in ab initio simulations is 
thus mandatory. On the atomistic scale, an explicit in¬ 
clusion of all solvent molecules in the simulation should 
be in principle the natural way to account for solvent 
effects. Due to the very large number of water and possi¬ 
bly other molecules required, this approach enormously 
increases the computational cost and limits at the same 
time the size of the system contained in the explicit di¬ 
electric medium^. The study of the solute-solvent inter¬ 
action at length scales larger than the molecular sizes 
would become virtually impossible. Investigations like 
structure predictions' or reaction path determination’ 
would become unaffordable in such a purely atomistic 
approach. Moreover, fully atomistic simulations of solva¬ 


tion effects would need to deal with the extensive sam¬ 
pling, required to characterize liquid configurations, and 
with the well-known limitations that current state-of-the- 
art ab initio methods present in describing liquid water, 
in particular regarding its structural and dielectric prop¬ 
erties. 

An implicit inclusion of the solute-solvent interactions 
could solve these issues. Starting from the earliest work 
of Onsager*’, the quantum chemistry community investi¬ 
gated implicit solvation models^'® extensively. In these 
models the solvent is introduced as a continuous homoge¬ 
neous and isotropic medium fully described by a dielec¬ 
tric function. Among them, the polarizable continuum 
model (PGM) developed by Tomasi and co-workers^>® is 
one of the most popular. In this approach a dielectric cav¬ 
ity surrounding the atomistic system is introduced where 
the permittivity takes on the value of one in regions oc¬ 
cupied by atoms and some different value characteristic 
for the dielectric solvent medium considered outside. 

Density functional theory (DFT) is a widely used 
method to investigate material properties at the atom¬ 
istic scale. In such ab initio calculations the electronic- 
structure problem is solved by minimizing the total en¬ 
ergy of the system which is a functional of the electronic 
density 

E[p\=T[p\+ j u(r)p(r)dr-I- ^ j p{Y)(j)[p\dY + 

( 1 ) 

where the four terms on the right side of Eq. (1) are, 
respectively, the standard kinetic energy, the interaction 
energy with an external potential, the electrostatic and 
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the exchange-correlation energy. For gas-phase molecu¬ 
lar simulations the potential ^(r) generated by a given 
charge density p(r) is given by the solution of the stan¬ 
dard Poisson (SPe) equation 

V^(/)(r) = —47rp(r). (2) 

An implicit inclusion of the solvent can be obtained by 
introducing a continuum dielectric cavity by means of a 
dielectric distribution e(r). Then the potential is given 
by the solution of the generalized Poisson equation (GPe) 


the non-linear Poisson-Boltzmann problem in a formula¬ 
tion which include ionic steric effects. The implemented 
algorithms take advantage of the chosen preconditioner 
for the minimization procedure. The accuracy of both 
GPe and PBe solvers has been tested for cases where 
an analytic solution is available. Two different methods 
have been implemented to describe the solvent surround¬ 
ing the atomistic system. We prove the effectiveness of 
our method in practical DFT calculations of electrostatic 
solvation energies of various test systems. 

II. GENERALIZED POISSON EQUATION 


V • e(r)V0(r) = —47rp(r). (3) 

If the system is surrounded by an ionic solution, an 
extra-term has to be added on the right side of Eq. (3). 
It accounts for the ionic distribution in the liquid and 
depends on the local electrostatic potential (/)(r). The 
resulting non-linear differential equation would therefore 
become 


V • e(r)V(/)(r) = -47r (p(r) -h (r)) , (4) 

where (r) is the local ionic concentration of the ions 
in the dielectric solvent, written as a sum of concentra¬ 
tion contributions Ci of ions of type i G {1, 2, and 

valence Zi, which in turn are (/>-dependent functionals: 

771 

[</'] (r) = gNa Z,Ci [(/)] (r), (5) 

i=l 

where e is the elementary charge and Na the Avogadro’s 
number. The most common expression for the Ci[(^] func¬ 
tional gives rise to the well-known Poisson-Boltzmann 
equation (PBe). 

Whereas several approaches and solvers exist for the 
gPeio-14, efficiejii; accurate solvers are still missing 
for the GPe and PBe cases. Some of these solvers handle 
sparse matrices obtained by low-order finite difference 
discretization of Eq. (3) , many of them keeping constant 
the permittivity in both inner and external regions of 
the dielectric cavity and then solving a standard Poisson 
problem. Furthermore, simpler and linearized forms of 
the Poisson-Boltzmann equation are usually considered, 
neglecting steric effects and overestimating ionic concen¬ 
trations close to highly charged surfaces and for multiva¬ 
lent ions. 

Fast and accurate GPe and PBe solvers which accu¬ 
rately works for a continuously varying dielectric func¬ 
tion could therefore play a key role in the extension 
of vacuum-based atomistic packages and allowing for 
quantum simulations in the presence of water, dissolved 
species, electrolytes, and non-aqueous solvents. 

In the present paper we present a minimization tech¬ 
nique which solves the generalized Poisson problem in 
some ten iterations of a SPe solver. In combination 
with a self-consistent procedure it enables us to solve 


As discussed above, the dielectric medium is described 
by means of a position-dependent dielectric distribution 
e(r). In order to solve numerically Eq. (3), both the 
potential </>(r), the charge density p(r) and the dielectric 
function e(r) are generally discretized on a finite grid. In 
principle also the generalized Poisson operator 


A = V-e(r)V 


( 6 ) 


should be discretized on the same mesh. It will be 
shown that depending on the adopted strategy to solve 
numerically Eq. (3), the discretization of the differen¬ 
tial operator A can be avoided in exchange of a iterative 
procedure based on a SPe solver. 

An alternative would be to solve the GPe iteratively as 
suggested in Ref. [ 15]. In this approach the polarization 
field introduced by the spatially varying dielectric func¬ 
tion is added as a source term to the charge density of 
the ordinary Poisson equation and the Poisson equation 
is solved repeatedly until self-consistency between the po¬ 
tential and the polarization charge density induced by it 
is reached. Our approach completes and simplifies this 
treatment, reducing considerably the number of SPe iter¬ 
ations, thereby presenting a robust and powerful iterative 
solver. 

Considering that reliable convergence can be an issue 
in mixing schemes, an alternative approach based on a 
minimization procedure is desireable. Such an approach 
can be based on an action integral whose Euler Lagrange 
equation is the GPe [Eq. (3)]: 


T = 


-V())(r)e(r)V^(r)) - 47rp(r)(^(r) 


dr . 


( 7 ) 


Any numerical minimization scheme can then be applied 
to solve the electrostatic problem. 

In Sec. IIB our strategy to solve Eq. (3), based on 
a preconditioned conjugate gradient (PCG) algorithm, 
will be presented. The preconditioner exactly represents 
the operator in the limit of a slowly varying dielectric 
constant and is based on the standard Poisson solver of 
BigDFT. 

The possibility of solving the generalized Poisson equa¬ 
tion under various boundary conditions (BC) will be very 
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Algorithm 1 Self-consistent iterative procedure (SC) 

1 Ater A 

1: Po =0 

2: for fc = 0,1, ■ • • do 

■ Pk =p/^ + pk 
4: solve 

5: Pl+i = i^Vlne- 

6 : Pl% = VPlAi + (1 - v)Pk^’~ 

rr „ Ater iter 

7: — Pfc+i — Pk 

8: end for 


important. In our approach the boundary conditions en¬ 
ter in a straightforward way by means of the precondi¬ 
tioner, i.e through the solution of the SPe. 

A. Self-consistent iterative procedure 


B. Preconditioned conjugate gradient 

Although a reasonably small number of iterations can 
be obtained in a self-consistent scheme, minimization 
techniques can produce more efficient methods to handle 
Eq.(3) if sophisticated schemes are utilized. In addition, 
the formulation as a minimization problem would allow 
to better control the convergence behavior. 

Solving this equation with a preconditioned steepest 
descent (PSD) method is essentially identical to the self- 
consistency approach of Sec. II A, once a standard Pois¬ 
son solver is taken as preconditioner. In particular, a 
good preconditioner for a PSD minimization, which in¬ 
verse applied to a residual vector provides the precon¬ 
ditioned residual Vk, is as follows: 

V^^Vkir) = e(r)V^z;fc(r) = -47rrfe(r). (11) 

Being the Laplacian of t'fe(r) related to the residual 
vector Vk by means of Eq. (11), the generalized Poisson 
operator becomes 


A strategy to solve the GPe for a given charge den¬ 
sity p(r) is by means of a self-consistent (SC) iterative 
procedure^^. Applying simple algebraic manipulations, 
Eq. (3) can be rewritten as 


V^(()(r) = —47r 


Xr) 

.eW 



=-47r [p(r)-f pP°'(r)] , 


( 8 ) 


where p^°^(r) is the polarization charge density. In this 
approach an extra-term p*‘®’’(r) 

p**-(r) = Xvine(r).V0(r) (9) 

induced by the spatially-varying dielectric function e(r) 
is added as a source to the charge density of the ordinary 
Poisson equation. Hence the GPe can be solved by a self- 
consistent loop on the potential (^(r), obtained by a SPe 
solver onto the second member of Eq. (8). The residual 
Tk , quantifying the convergence, is the difference between 
the extra terms of Eq. (9) between subsequent iterations. 
Algorithm 1 describes the procedure. In order to stabilize 
the iterative method, a linear mixing of the extra-term 
^zter(r) steps fcth and {k l)th has been introduced 
tuned by means of the mixing parameter r]. 

The polarization charge induced in the dielectric 
medium can be easily related to the extra-term of Eq. 
(9): 

/°'(r)=X-(r) + "^^p(r). (10) 

e(r) 

This charge represents the response of the surrounding 
implicit dielectric. It lies in the transition region between 
the inner and outside part of the cavity and stabilizes the 
solute density enveloped by the solvent. 


Avkir) = \/■ e{r)\/vk{r) 

= Ve(r) • Vnfe(r) - 47rrfc(r). 

Such a PSD approach can be described by Algorithm 
2 with Pk = 0. Fixing = 1 it corresponds to the 
previously described self-consistent approach. 

In a PSD scheme the number of iterations I needed 
for convergence is proportional to the condition number 
K (i.e. the ratio between the largest and the smallest 
eigenvalue of the product operator V~^A). Minimiza¬ 
tion methods with a faster convergence rate than the pre¬ 
conditioned steepest descent algorithm can significantly 
improve the convergence speed. 

We use a preconditioned conjugate gradient scheme, 
where I cx -A. In such minimization procedure a good 
preconditioner can lower k and, therefore, the overall 
number of iterations. Algorithm 2 describes the imple¬ 
mented PCG procedure to compute the electrostatic po¬ 
tential (/)(r) starting from a given charge density p{r). 
The minimization procedure starts from an initial gradi¬ 
ent To computed on an input guess po. V is the precon¬ 
ditioner which inverse has to be applied to the residual 
vector Tfc returning the preconditioned residual and, 
finally, pk is the solution of Eq. (3). The convergence cri¬ 
terion is imposed on the Euclidean norm of the residual 
vector r/j. 

Both the performance and accuracy in a PGG scheme 
critically depend on the preconditioner chosen. We im¬ 
plemented a preconditioner based on the solution of the 
standard Poisson equation (namely on a standard Posson 
solver). Once a residual vector is given at the step 3 
of Algorithm 2, we define a preconditioned residual from 
the following equation: 

V^^Vkir) = vXOV^X(r)\/e(iO] = -47rrfc(r). (13) 
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Algorithm 2 Preconditioned conjugate gradient 
ro = —47rp — A(j)o, P-i = 0 


for fc = 0,1, • • • do 

Vk = 

Pk =Vk+ PkPk-1 (where h = , fc 7 ^ 0 ) 

"'= = (k%i) 

(('fe+l =((>*;+ OfePfe 
Tfc+i =rk- akApk 


8: end for 


This equation has to be solved with respect Vk (r) once 
rfc(r) is given. 

In addition to speeding up the PCG procedure, the 
preconditioner defined by Eq. (13) retains a further fea¬ 
ture which guarantees accuracy and fast performance for 
the whole electrostatic solver. In step 7 of Algorithm 
2 we have to apply the generalized Poisson operator A 
to the preconditioned residue pk, which means, thanks 
to step 3, app lying it to Vk- Using a change of variable 
v'j.{r) = ^^/e{r)vk{Y), the GPe becomes 

V • e(r)Vz;fc(r) = \A(r)V\fe(r) - Vk{r)\/^y/^. (14) 

Now simple algebraic manipulations and Eq. (13) al¬ 
low to rewrite the generalized Poisson operator A as 

Avkir) = ■ e{r)\/vk{r) (15) 

= -Vk{r)q{r) - 47rrfe(r), (16) 

where q{r) = ^e(r)V^-\/e(r) is calculated once at the 
beginning of the PGG procedure and kept fixed for the 
whole minimization loop. Therefore thanks to the cho¬ 
sen preconditioner, the action of the operator A can be 
simplified to a simple multiplication between the po¬ 
tential ffe(r) and a constant vector q(r) related to the 
spatially-varying dielectric function e(r). This feature, 
which provides the exact operator output, makes our 
PGG procedure robust and fast, avoiding any finite dif¬ 
ference differentiation. Eurthermore, reducing the PGG 
algorithm to simple vector operations, makes its paral¬ 
lelization straightforward, delegating it to the chosen SPe 
solver. A similar discussion holds for the boundary con¬ 
ditions, which enter in a natural way by means of the 
preconditioner, i.e through the solution of the ordinary 
Poisson equation, both in the SG and PGG algorithm. 



c 

_o 

o 


U 

Q 


FIG. 1. Analytical three dimensional functions used as bench¬ 
mark fields for both SC and PCG solvers along a particular 
direction passing though the box center and parallel to the 
y axis. Red dash line: potential (/)(r); red dot line: charge 
density p(r); black solid line: dielectric function e(r). 


Algorithm 2) , we used the Interpolating Scaling Function 
(ISF) Poisson Solver, allowing to obtain highly accurate 
electrostatic potentials for free, wire, surface, and peri¬ 
odic boundary conditions at the cost of 0(Nlog(N)) op¬ 
erations, where N is the number of discretization points 
(see Ref. [ 14]). 

To test both solvers analytic three dimensional func¬ 
tions have been used. An orthorhombic grid of uniform 
mesh spacing hg^id and (nx,ny,nz) points in each direc¬ 
tions has been used. Fig. 1 shows plots of these bench¬ 
mark fields along a particular direction passing through 
the box center and parallel to the y axis. All functions 
depend on the radial distance r from the center of the 
simulation domain. 

A normalized Gaussian function has been chosen for 
the electrostatic potential (j){r) (red dash line in Fig. 1) 
and the charge density p(r) has been derived from the 
chosen potential and dielectric function, applying the 
generalized Poisson differential operator of Eq. (6) (red 
dot line). In order to reproduce the dielectric environ¬ 
ment in electrostatic problems typically used for a solute 
system embedded in a solvent (i.e. a cavity where the ma¬ 
jority of the atomic charge density is confined), the error 
function 1 -|- (cq — l)h{do,/S.]r) has been chosen to rep¬ 
resent the spatially varying electric constant e(r) (solid 
black line in Fig. 1), where 


C. Numerical results 


h{do,A;r) 


1 

2 


1 -b erf 


r-do 

A 


(17) 


Both the self-consistent iterative procedure (Algorithm 
1) and the preconditioned conjugate gradient minimiza¬ 
tion scheme (Algorithm 2) have been implemented and 
tested. As SPe solver (step 4 of Algorithms 1 and step 3 of 


Here A is a parameter which controls the transition re¬ 
gion (« 4A wide) between the inner part within the cav¬ 
ity of radius do and the external part, and Cq is the dielec¬ 
tric constant of the surrounding medium. These bench- 
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Free Surface Periodic 



FIG. 2. Euclidean norm of the residual vector rt (top graphs) and accuracy of the numerical solution (lower graphs) for SC 
(blue circles) and PCG (black squares) solvers with free, surface and period boundary conditions. 


mark functions have been used for both free, surface and 
periodic boundary conditions. 

To implement the nabla differential operator needed 
for Algorithm 1, central, forward and backward finite 
difference filters of order 16 have been used, which match 
the accuracy of the underlying SPe solver. We remark 
that the use of these finite difference filters is not needed 
for Algorithm 2 as soon as the correction term of Eq. (16) 
is pre-calculated. 

Fig. 2 shows solver performances. The top graphs re¬ 
port the residual norm as function of the iteration num¬ 
ber. The residual norm is defined has the Euclidean norm 
of the residual vector r^. Graphs in lower panel present 
the output accuracy as a function of the iteration num¬ 
ber. The accuracy in the whole paper is dehned as the 
maximum value of the difference between the final nu¬ 
merical solution and the analytical potential. In this test 
case a cubic box of length 10 [arb. units] has been chosen 
with Ux = Tiy = Uz — 300. The Gaussian variance for 
the potential ^(r) was a = 0.5, and the parameters of 
e(r) were do = 1-7, A = 0.3 and eg = 78.36 (all in [arb. 
units]). The mixing parameter in step 6 of Algorithm 


1 has been hxed to be ry = 0.6, resulting in a robust 
convergence for all cases. Lower values slow down the 
convergence for the chosen test functions. 

The PGG solver (black squares) exhibits a faster con¬ 
vergence with respect to the SG one (blue circles), reach¬ 
ing an accuracy of ^ 10“^° with some ten iterations. Fur¬ 
thermore its behavior does not change with the bound¬ 
ary conditions as is the case with the SG algorithm. It 
is worth remembering that each PGG iteration involves 
only a single solution of the ordinary Poisson equation 
and as well as of fully parallelizable vector operations. If 
an accuracy of ^ 10“"^ — 10“® is enough, then some five 
iterations solve the electrostatic problem. These features 
make the developed PGG algorithm together with the 
chosen preconditioner of Eqs. (13) very efficient for atom¬ 
istic calculations where the generalized Poisson equation 
needs to be solved repeatedly. Performances of the im¬ 
plemented PCG procedure are also higher than multigrid 
approaches to solve the GPe, where a number of itera¬ 
tions between 17 and 25 are needed to reach an accuracy 
of ~ 10-®i'L 

For the sake of completeness, the preconditioned steep- 
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est descent scheme (Algorithm 2 with jSk = 0) has been 
tested with the preconditioner described by Eq. (11). 
Residual norm convergence as well as the accuracy of the 
solution (not reported) behaves like the self-consistent 
approach (Algorithm 1, blue circles of Fig. 2). An inte¬ 
gration of DIIS in the PSD loop effectively lowered the 
iteration numbers, but didn’t show better performances 
with respect to the PCG approach. 


III. POISSON-BOLTZMANN EQUATION 


The generalized Poisson equation so far discussed holds 
for a solute plunged in a neutral solution where no mobile 
charges are present. In order to extend the library to 
ionic solutions, effects of mobile ions have to be taken into 
account which, being free to move inside the dielectric 
medium, modify the spatial distribution of charge and 
potential close to the interface and giving rise to the well 
known double layer. 

In general, mobile charges can be included in the elec¬ 
trostatic problem by means of a continuum mean-field 
approach, assuming point-like ions in thermodynamic 
equilibrium. Once the equilibrium is reached, ionic con¬ 
centrations explicitly depend on the local electrostatic 
potential (f>{r). Following this idea, the potential ^(r) 
generated by a charge density p(r) placed in contact with 
a ionic solution can be extracted solving the generalized 
Poisson equation (Eq. (4)). Several models have been 
proposed in the literature for the ionic bulk concentra¬ 
tions [see Eq. (5)]. Gouy^' and Ghapman^® proposed a 
Boltzmann distribution 


Ci[(t)]{r) 


= Ci exp 


\ kT 


(18) 


where k is the Boltzmann constant and T the absolute 
temperature of the solvent. Gombining Eqs. (4,5,18) the 
well-known Poisson-Boltzmann equation can be recov¬ 
ered. It arises from the equilibrium between thermal and 
electric forces, which depend, respectively, on the ionic 
concentration and electrostatic potential. At equilibrium 
the total average force must be zero and Eq. (18) holds. 
In the regions where the electrostatic energy is smaller 
then kT, i.e. Zie(j){r)/kT <C 1, the exponential of Eq. 
(18) can be approximated with a linear function of (()(r), 
switching the non-linear problem of Eq. (4) into a linear 
one. 

The Poisson-Boltzmann equation correctly predicts 
ionic proHles close at the solid-liquid interface with ionic 
solutions. However it strongly overestimates ionic con¬ 
centrations close to highly charged surfaces or multiva¬ 
lent ions. In order to overcome these drawbacks, several 
models have been proposed^®. Finite ion size effects can 
be included in the model by introducing an additional in¬ 
ternal force. Using a Bikerman-type expression to model 
steric effects^'^ and imposing that at equilibrium the to¬ 
tal average force (thermal, electric and steric) must go to 


zero, a Langmuir-type equation for the ionic concentra¬ 
tions can be found 


Ci[(t)]{r) 


c^exp 


Zie(j)(v)\ 
kT ) 


1 + 


— 

E S' 

„max 
3 = 1 J 


exp 


kT ) 



In Eq. (19) are the maximum local concentration 
that a ionic species of effective radius Ri can attain. This 
last is given by 


C = -, (20) 

-ttRINa 

where p is the packing coefficient. The combination of 
Eqs. (4,5,19) give rise to the so-called modified Poisson- 
Boltzmann (MPBe) equation^It accounts for hnite 
ion size effects representing an extension of the PBe, pre¬ 
venting ion concentrations to exceed Note that if 

^max QQ Yj g {l,2,...,m}, Eq. (19) is reduced to 
Eq. (18) and the standard Poisson-Boltzmann equation 
holds. 

Both PBe and MPBe equations consider ions in so¬ 
lution as pointlike. A further extension can account 
for finite size effects in an explicit way, describing ions 
by means of insulating spheres plunged in the dielectric 
solvent^”*. This correction allows for a local modihca- 
tion of the solvent permittivity as well as the inclusion of 
two new forces: one related to a non-uniform dielectric 
which tend to move ions into high-permittivity regions; 
another due to the interaction of the ion dipole and the 
non-uniform local electric field (dielectrophoretic force). 

The Poisson-Boltzmann equation is more difficult to 
solve than the generalized Poisson equation due to 
the fact that it leads to a non-linear problem. If 
Zie(j){r)/kT I, Eq. (18) can be approximated by a 
linear function of (j){r) , and the Poisson-Boltzmann prob¬ 
lem of Eq. (4) becomes linear. In Sec. IIB a particular 
preconditioner, i.e. Eq.(13), has been proposed for the 
solution of the GPe. Furthermore the generalized Pois¬ 
son operator has a linear term with respect to (j){r) [Eq. 
(16)]. Therefore Algorithm 2 is expected to solve the lin¬ 
ear regime of the PBe, where the A operator is now given 
by Eq. (16) with an additional linear term: 

Avk{r) = V • e(r)Vufc(r) -k47rp“”®[(/)](r) 

= -Vkir) ^g(r) -b ^ - AirrAr) , 

( 21 ) 

where q{r) has been defined in Sec. IIB. In the gen¬ 
eral non-linear case, at variance to the generalized Pois¬ 
son equation discussed in the previous sections, there 
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Algorithm 3 Self-Consistent iterative procedure for 
the Poisson-Boltzmann equation 

1: set Po 

2: for fc = 0,1,... do 

o „tot I Aons 

3: Pk = P + Pk 

4: solve V • eV(j>k = —(Algorithm 1 or 2) 

5: compute pl°A = [<(>&] 

: Pk+i = VPk+1 + (1 - V)Pk 

PT „ Aons Aons 

7: Tk+i = Pk+i - Pk 

8: end for 


is no general consensus on what flavor of the Poisson- 
Boltzmann equation is most appropriate for various prob¬ 
lems. Therefore a (M)PBe solver should allow to solve 
various formulations of this equation. Since it is unlikely 
that for all variants an action integral can be established 
(allowing for a minimization scheme), a self-consistent 
approach should be the most appropriate algorithm. 

Such a procedure has been implemented and it is de¬ 
tailed in Algorithm 3. The ionic charge density 
is included as source term to the charge density of the 
generalized Poisson equation which, in turn, is solved re¬ 
peatedly until self-consistency between the potential and 
the ionic charge induced by it is reached. Starting with 
an initial input guess for the ionic charge density, 

on each self-consistent cycle a generalized Poisson solver 
(Algorithm 1 or 2) is applied to numerically compute 
the electrostatic potential (/)(r). Then the ionic charge is 
computed using the new potential and mixed with the 
old density by means of a linear scheme tuned by the 
parameter r]. 

In order to speed up performances of the (M)PBe 
solver, an improved version is reported in Algorithm 4. 
On each iteration k the electrostatic problem at step 4 
is solved only for the previous residual vector once 
pions been updated at step 3. Then the overall solu¬ 
tion is given at step 5 as sum over all potential corrections 
(p'f,. This procedure substantially corresponds to the gen¬ 
eral self-consistent approach of Algorithm 3, now using 
on each step information of the previous as input guess. 
It is worth noting that the improved Algorithm 4 can be 
coupled only with the faster PCG solver (Algorithm 2) 
for the GPe. 


A. Numerical performances 

To show performances and accuracy of the Poisson- 
Bolzmann solver, both Algorithms 3 and 4 have been 
applied to analytical test cases. Following the strategy 
reported in Sec. IIC, all the involved helds and functions, 
i.e. the potential, the charge density and the operator 


Algorithm 4 Improved Self-Consistent iterative 
procedure for the Poisson-Boltzmann equation 

1 ^^4- r\ Aons r\ Aons A 

1: set (pi — 0, Po — 0; Pi 5 ^0 ^ P 

2: for k = 1,... do 

3 Aot A I Aons Aons 

■■ Pk =r^_i+Pk - Pk-i 
4: solve V -tVcfif. = (Algorithm 2 with a residual 

vector r(.) 

5: (j^k = 4^k (j^k 

6: compute 

7: Pl°+T = VPk+l + (1 - V)pr‘ 

8 „ Aons Aons 

■ rk+i = Pk+i - Pk 

9: end for 


have been discretized on the orthorhombic three dimen¬ 
sional grid and same parameters have been used to set up 
all analytical functions. The electrostatic problem lies in 
a cavity where the great majority of the charge density is 
confined, described by means of a dielectric error function 
e(r) (solid black line of Fig. 1). A normalized Gaussian 
function has been chosen for the electrostatic potential 
(/)(r) (red dash line of Fig. 1), whereas the charge den¬ 
sity p(r) has been derived from the chosen potential and 
dielectric function, applying the Poisson-Boltzmann dif¬ 
ferential operator. 

A monovalent {Zi = —Z 2 = 1) binary aqueous elec¬ 
trolyte solution has been considered, with a close packing 
coefficient p = 0.74, effective ionic radius = 3 [A] and 
bulk ion concentrations = 100 [mol/m^] kept fixed for 
all ions. As discussed in Sec. Ill, the PGG algorithm has 
been applied to solve the linear Poisson-Boltzmann equa¬ 
tion with the same GPe preconditioner of Eq. (13). In 
Fig. 3 the euclidean norm of the residual vector (black 
squares) and the accuracy of the numerical solution (blue 
circles) have been reported. Similar performances have 
been found with respect to the generalized Poisson solver, 
reaching an accuracy of ^ 10“^'^ with some ten iterations 
and proving that the PCG procedure is a well suited and 
fast method also for the linear regime of the PBe. 

In the general non-linear case, Algorithms 3 and 4 have 
been applied. An input guess Pq°'^^ = 0 has been chosen. 
We found that only a relatively small number of self- 
consistency iterations is needed in this approach to solve 
both the Poisson-Boltzmann and the modified Poisson- 
Boltzmann problem. Fig. 4 shows the Euclidean norm of 
the residual vector (black squares) and the accuracy 
of the numerical solution (blue circles) for the modified 
Poisson-Boltzmann solver with free boundary conditions. 
The solver exhibits a fast convergence, reaching an accu¬ 
racy of ~ 10“^° with some hve iterations. The inset 
shows as given by Eqs. (5,19), revealing an ionic 

density saturation in the solvent for electrostatic poten- 
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FIG. 3. Euclidean norm of the residual vector (black 
squares) and accuracy of the numerical solution (blue circles) 
for the linear Poisson-Boltzmann equation with free boundary 
conditions (Algorithm 2 has been applied with A the linear 
Poisson-Boltzmann operator). 



FIG. 4. Euclidean norm of the residual vector Vk (black 
squares) and accuracy of the numerical solution (blue circles) 
for the modified Poisson-Boltzmann solver with free bound¬ 
ary conditions. Numbers between round brackets represent 
the number of iterations needed to solve the GPe at step 4 of 
Algorithm 4. The inset shows given by Eqs. (5,19). 


tials which absolute value is higher than ^ 0.01 [a.u.]. 
Numbers between round brackets represent the number 
of iterations needed to solve the GPe at step 4 of Algo¬ 
rithm 4 using the PCG scheme. The convergence crite¬ 
rion for the GPe solver has been fixed equal to the one 
of Algorithm 4, except for the first two iterations where 
very accurate potentials do not change the overall per¬ 
formances of the PBe solver. Furthermore its behavior 
does not change with the boundary conditions which are 
managed by means of the generalized Poisson solver. The 
latter also fixes the final accuracy of the self-consistent 
procedure. It is worth noting that performances and par¬ 
allel efhciency as well as boundary conditions of both GPe 
and (M)PBe solvers are, eventually, delegated to the un¬ 
derlying SPe solver. 

IV. ELECTRONIC STRUCTURE 
COMPUTATIONS 

Effects of complex wet environments surrounding an 
atomistic system can be approximately included into den¬ 
sity functional calculations by simply introducing a di¬ 
electric cavity surrounding the atomistic system and tak¬ 
ing in Eq. (1) the electrostatic potential (j)[p[Y)\ as so¬ 
lution of the generalized Poisson or Poisson-Boltzmann 
equation. The PCG solver for the GPe (Algorithm 2) 
has been implemented in the electronic-structure package 
BigDFT^®, extending the capability of the code beyond 
vacuum-simulations. Two distinct approaches have been 
implemented and tested to build up a dielectric cavity 
enveloping the atomic system. In both approaches the 
cavity, mapped by the dielectric function e(r), is fully 
differentiable and continuous in the whole simulation do¬ 


main. 

In the first approach a function e(r) is defined starting 
from spherical-symmetric atom-centered cavities. Each 
sphere depends only on the radial distance with respect 
to a fixed atomic position. The whole rigid cavity is kept 
fixed during the whole SCE cycle in a DFT simulation. 

On the other hand, it could be argued that regions oc¬ 
cupied by atoms and their associated electronic charge 
density are strictly related. In other words, the actual 
value of the charge density might determine how much 
“space” is occupied by the solute. Starting from this idea, 
a cavity can be build up directly from the electronic den¬ 
sity, as it has been shown in various publications^^’^®’^^. 
In this second approach the dielectric function is not ex¬ 
plicitly space-dependent, but can be implicitly mapped 
by means of the electronic charge density. 


A. Rigid cavity 

The most widespread continuum solvation model, 
which tries to include the effects of a surrounding dielec¬ 
tric medium in an implicit way, is the polarizable contin¬ 
uum model developed by J. Tomasi and co-workers^’*’^*. 
In the PCM formulation the cavity surrounding the so¬ 
lute is sharp and discontinuous, and a polarization charge 
density is exactly localized at the vacuum-dielectric in¬ 
terface. In this way the dielectric environment is repre¬ 
sented by an effective surface polarization charge, reduc¬ 
ing the electrostatic problem into a two-dimensional one. 
Furthermore the cavity can be considered rigid since it 
depends only on the atomic coordinates which does not 
vary during a SCF cycle in DFT simulations. 

For a first test of the electrostatic solvers (Algorithm 
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1 and 2) integration in an electronic-structure code, 
a PCM-like cavity has been considered. The dielec¬ 
tric function e(r) is given by the product of spherical- 
symmetric atom-centered error functions. Differently 
from the early PCM model, we here define a cavity which 
is fully differentiable and continuous. In particular, for a 
system of N atoms of coordinates R, (for i = 1,..., N) 


e(r, {RJ) = (eo - 1) |]J A; d{r, R*)) | (^2) 

where eo is the dielectric constant of the surrounding 
medium and the function h is defined by Eq. (17). In Eq. 
(22) d(r,R,) = ||r-R,||, d i are the cavity radii which 
depend on the particular atom species, and A a parame¬ 
ter (kept fix for all atoms) which controls the transition 
region (« 4A wide) from 0 to 1 where the polarization 
charge is located. Starting from Eq. (22), all vectors 
which explicitly depend on e(r) can be analytically com¬ 
puted. The cavity is uniquely defined once Ri, di and 
A are fixed. All environment-dependent fields are cal¬ 
culated once at the start of the solver and kept fixed 
throughout the procedure. 

In must be noticed that this definition of the cavity 
relies on the explicit dependence of e(r, {Ri}) from the 
atomic coordinates R^ (and, consequently, of the system 
Hamiltonian). This dependence gives rise to additional 
terms when atomic forces are computed. The analytical 
rigid cavity above described should overcome this prob¬ 
lem allowing a direct analytic calculation of such addi¬ 
tional contributions. Moreover, the values of di and of 
A have to be tuned by the user, usually by choosing a 
solvent-dependent scaling factor with respect to empiri¬ 
cal Van der Waals radii^'"*. A solution that would remove 
part of this arbitrariness should therefore avoid an ex¬ 
plicit use of atomic coordinates in the cavity mapping. 
This will be the subject of the following section. 


Several features make the self-consistent cavity more 
advantageous with respect to the rigid one. First, once 
the electron charge density is given, only two parame¬ 
ters uniquely define both the cavity and the transition 
region for the whole atomic system. Furthermore, since 
the dielectric function does not explicitly depend on the 
atomic coordinates, the evaluation of ionic forces can be 
done without modifications with respect to a simulation 
in gas-phase. 

Among several ways to define the functional depen¬ 
dence on the electronic charge density, the self-consistent 
continuum solvation (sees) model developed in Ref. [ 
15] has been implemented. It allows to fit experimental 
solvation energies on a set of 240 neutral solutes with a 
mean absolute error of 1.2 kcal/mol: 



1 

„ele.c „ 

H f-'max 


e(p"'“) = < 


Pmin ^ P ^ ^ Pmax ’ 

(24) 


eo 

V. 

„ele.c ^ „ 

H ^ h'min 


where w{x) is a continuous smooth function describing 
the transition region between vacuum (where atoms are 
placed) and the full dielectric medium; 


(x) - sin(z(a:))] ; 

(25) 

In 

zix) = 2 tt — 


(26) 


A good description of this region by means of Eq. (24) 
is mandatory for the procedure convergence as well as 
for the mapping of the polarization charge there con¬ 
fined. Since e(r) explicitly depends on p®*®'^(r), its vari¬ 
ation needs to be included in the Kohn-Sham (KS) po¬ 
tential. Starting from Eq. (3) and integrating by parts, 
the electrostatic energy can be rewritten as 


B. Charge-dependent cavity 

For this definition of the cavity, the dielectric function 
does not explicitly depend in the atomic positions, but 
implicitly via the charge density 


Ees[p] = ^J p(/>[p]dr = ^ J e[p](V((i[p])^dr. (27) 

Its functional derivative with respect to p gives the 
electrostatic potential (j) plus an additional term r'e(r) 


e(r)=e[p^'-](r). (23) 

This approach allows the cavity to self-consistently 
change during the SCF loop, strictly following the mod¬ 
ification of the electronic charge density. A cavity sur¬ 
rounding an atomic system can be generated by means 
of two threshold parameter, i.e. Pmax and Pmin, fixing 
e(r) = 1 in regions when p®*®‘^(r) > pmax and e(r) = eo 
when p®'®^(r) < Like dt in the rigid case, pmax 

fix the width of the cavity, whilst the extension of the 
transition region, previously defined by A, is now tuned 
Ly Pmin- 


Ve{r) 


1 de(p®'“(r)) 
Stt 


|V(/.(r)|2 , 


which has to be added to the KS potential. 


(28) 


C. Solvation Free energies 

In order to test the integration and performances of the 
generalized Poisson solver into ab initio codes, the whole 
procedure previously described, i.e. Algorithm 1 and 2 of 
Sec. IIA and IIB, rigid and charge-dependent cavities of 
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Sec. IV A and IV B as well as the additional term of Eq. 
(28) have been integrated in the main BigDFT package^^. 
This extension allows to handle complex wet environ¬ 
ments in electronic-structure calculations, including im¬ 
plicitly the effects of a solvent surrounding an atomic 
system. 

The electrostatic solvation energy is defined as differ¬ 
ence between the total energy of a given atomic system in 
the presence of the dielectric environment and the energy 
of the same system in vacuum 

A(^e/ ^ Qel _ qO ( 29 ) 

A full comparison with experimental solvation energies 
needs the inclusion of non-electrostatic contributions. In 
this case the main terms in the solute Hamiltonian, as 
introduced by PCM®, are 

-f + G^*" -k G*"* -k PAV, (30) 

where AG®* is the electrostatic contribution, G®“'“ the 
cavitation energy, i.e. the energy necessary to build up 
the solute cavity inside the solvent medium. G®®^ is a re¬ 
pulsion term representing the continuum counterpart of 
the short-range interactions induced by the Pauli exclu¬ 
sion principle, whilst G'**® reflects van der Waals inter¬ 
actions. The thermal term G*™ accounts for the vibra¬ 
tional and rotational changes and, finally, PAV includes 
volume change in the solute Hamiltonian. 

The inclusion of all non-electrostatic contributions 
goes beyond the aim of the present paper, where test¬ 
ing of the generalized Poisson solver and its integration 
into first principle atomistic calculations is the main goal. 
Therefore only electrostatic solvation energies have been 
computed for a set of small neutral organic molecules. 

Water has been chosen as solvent for all cases, with a 
dielectric constant of eg = 78.36 (experimental value at 
low frequency and ambient conditions). Final energies for 
all molecules have been extracted after a full geometry 
optimization both in vacuum and in aqueous solution. In 
all cases the surrounding dielectric medium lowers the to¬ 
tal energy of the system with respect to vacuum, because 
the polarization of the dielectric stabilizes the solutes. 

Table I reports electrostatic solvation energy AG®* ob¬ 
tained both with rigid and self-consistent cavities under 
free boundary conditions. As previously stated, a criti¬ 
cal point of PCM approaches is the choice of shape and 
size of the cavity, which should mimic the solute incor¬ 
porating the whole atomic charge density. From its first 
formulation^®, PCM atomic radii Ri were fixed propor¬ 
tional to the van der Waals radii, Ri = f ■ In 

our rigid model, a proportional factor of / = 1.2 has 
been fixed and the Pauling’s set of atomic radii has been 
considered®**"®^ (except for the hydrogen atoms bound to 
heteroatoms which have a radius value of 1 A). Having 
a further degree of freedom with respect to sharp PCM 
cavities, a A = 0.5 [a.u.] has been tuned and kept con¬ 
stant for all atoms. 


TABLE I. Electrostatic solvation energies AG®* (in kcal/mol) 



PCM’^ 

rigidBigDFT 

SCCSqe® 

SCCSbikOFT 

NHg 

-6.65 

-6.28 

-5.39 

-5.35 

HaO 

-8.98 

-8.32 

- 8.21 

-8.23 

CH 4 

-0.61 

- 1.20 

- 0.68 

-0.63 

CH 3 OH 

-6.78 

-6.57 

-5.89 

-5.83 

CH 3 NH 2 

-4.51 

-5.71 

-4.53 

-4.45 

CH 3 CONH 2 

-12.53 

-12.97 

-11.87 

-11.87 


^ Polarizable continuum model results obtained with GAUSSIAN 
09 ’'^ and Pauling’s set of atomic radii'^^ 

^ Self-consistent continuum solvation model results from periodic 
BC calculations in Ref.^ ’ 


As reference, in Table I sees calculations from Ref. [ 15] 
together with PCM calculations have been reported. The 
latter were performed with the GAUSSIAN 09 code®®, but 
using the version of PCM which was the default in GAUS¬ 
SIAN 03, as specified by the keyword g03defaults. Only 
electrostatic solvation effects were included in the calcu¬ 
lation and, to simplify the comparison, the simple van der 
Waals surface was adopted with Pauling’s set of atomic 
radii®* (explicit hydrogens), but without the additional 
smoothing used to describe the solvent-excluded surface. 
The Perdew-Berke-Ernzerhof (PBE) functional ' and the 
extended triple-zeta 6-311-|-g(d,p) basis set were used for 
both geometry optimizations and energy calculations, in 
vacuum and in solution, consistently with the set up used 
for the parameterization of the electrostatic solvation en¬ 
ergy in sees*®. 

In all calculations soft norm-conserving pseudopoten¬ 
tials including non-linear core correction®®’®** along with 
PBE functional were used to describe the core electrons 
and exchange-correlation, respectively. 

For the self-consistent cavity, values of Pmax = 5 • 10"® 
and Pmin = 1 • 10"® have been fixed which produces a 
mean absolute error of 1.20 kcal/mol for the solvation 
energies of a database of 240 molecules*®. A perfect 
agreement has been reached which confirms the perfor¬ 
mance and reliability of the integrated generalized Pois¬ 
son solver. 

In order to test and validate the whole library in 
BigDFT, the effects of grid resolution and boundary con¬ 
ditions on the electrostatic solvation energies for the same 
set of neutral molecules have been investigated. The 
global accuracy is strictly related to the size of the sim¬ 
ulation box and the spatial grid resolution ft-grid- How¬ 
ever a decrease of the latter can affect the whole cost of 
the calculation both in time and memory usage. Conse¬ 
quently it is worth to investigate the effects of the dielec¬ 
tric medium’s inclusion in a DFT run with respect to the 
vacuum case. 

In this respect Kohn-Sham total energies have been ex¬ 
tracted both in vacuum and in the presence of the solvent 
(water) as function of the spatial grid hg^id. Its accuracy 
is reported in Fig. 5 as difference with the reference value 
at /igrid = 0.20 bohr. Results show that the accuracy is 
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not affected by the surrounding dielectric environment 
(filled symbols, solid lines) with respect to the vacuum 
case (empty symbols, dash lines). Following the same 
guidelines, the accuracy of the electrostatic solvation en¬ 
ergy [Eq. (29)] with respect to the spatial grid /igod has 
been investigated and reported in Fig. 6 . Results show 
that a /igrid = 0.30 bohr in BigDFT with free boundary 
condition provides electrostatic solvation energies with 
an accuracy lower than ~ 10 “^ kcal/mol. 

Since isolated systems embedded in a wet environment 
are the subject of interest, spurious long range electro¬ 
static interactions with periodic images due to artificially 
imposed periodicities along certain direction can be prob¬ 
lematic. In the developed procedure the boundary con¬ 
ditions enter through the preconditioner, i.e. the SPe 
solver. Since in the BigDFT Poisson solver all the com¬ 
mon boundary conditions such as free, wire and surface 
are exactly implemented, such spurious interactions do 
not exist at any stage of our approach. Fig. 7 reports 
the difference between the electrostatic solvation ener¬ 
gies computed with periodic and free boundary condi¬ 
tions as function of the periodic cell length and keeping 
fixed the spatial grid hg^id = 0.20 bohr. Molecules have 
been ordered according to the strength of their electro¬ 
static dipole, from largest to smallest. Dipole values for 
each molecule are reported on Table II both in vacuum 
and in the present of a water solvent. It can be no¬ 
ticed that the presence of the polar solvent increases the 
electrostatic dipoles of the molecules^®. As it might be 
expected, the free BC calculation represents the asymp- 
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FIG. 6. Accuracy of the electrostatic solvation energy with 
respect to the spatial grid hgrid with free boundary conditions. 


totic results for periodic BC calculation of increasing box 
sizes. Interactions with image-molecules are less relevant 
for molecules with small dipole moment like CH 4 , but 
are not negligible for molecules with larger dipoles such 
as CH 3 CONH 2 . In such cases unrealistic large periodic 
boxes are required for periodic BC to reproduce the free 
boundary condition results with high accuracy. 

Numerous processes of practical interest involve sur¬ 
faces in contact with neutral or ionic solvents, leading to 
an induced polarization charge of the dielectric medium 
or an electric double layer. The BigDFT package allows 
to use exact surface boundary conditions avoiding spuri- 



FIG. 5. Accuracy of the Kohn-Sham total energy with re¬ 
spect to the spatial grid /igrid- Galculations for all molecules 
have been done both in vacuum (empty symbols, dash lines) 
and with a surrounding dielectric environment (filled symbols, 
solid lines) in free boundary conditions. 


FIG. 7. Difference between the electrostatic solvation energy 
computed with periodic and free boundary conditions as func¬ 
tion of the periodic cell length. Molecules have been ordered 
as function of their electrostatic dipole norm, from largest to 
smallest (see Table II). 
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TABLE II. Molecule dipole norm in vacuum and water solvent 
(in Debye). 



Dipole’'“‘=““™ 

Dipole”“*=” 

CH 3 CONH 2 

3.88 

5.76 

H 20 

1.81 

2.41 

CH 30 H 

1.57 

2.14 

NH 3 

1.49 

1.98 

CH 3 NH 2 

1.27 

1.78 

CH 4 

0.00 

0.00 



FIG. 8 . Ti 02 surface in contact with water: isosurfaces of 
the polarization charge in the implicit dielectric medium. 


ous interaction in the direction orthogonal to the surface 
exposed to the wet environments. To show a further ap¬ 
plication of the solvation library in BigDFT, a Ti02 sur¬ 
face in contact with pure water has been simulated. The 
full DFT simulation in presence of the solvent has been 
initialized starting from its relaxed state in vacuum. Fig. 
8 shows the Ti 02 wet surface as well as isosurfaces of the 
polarization charge density as given by Eq. (10). In 
this test case, we found that the contact with the dielec¬ 


tric medium provides an electrostatic solvation energy of 
-46.60 kcal/mol. 

V. CONCLUSIONS 

In the present work a library to handle complex 
wet-environments in electronic-structure calculations has 
been presented. It allows to include on the atomistic scale 
effects of an aqueous environment in an implicit way. 
The solver is able to handle both the generalized Poisson 
and several variants of the Poisson-Boltzmann equation. 
The core of the generalized Poisson solver is a precon¬ 
ditioned conjugate gradient algorithm which allows to 
numerically solve the minimization problem with some 
ten iterations. The same algorithm works for the linear 
case of the Poisson-Boltzmann equation, whilst for the 
general case a self-consistent procedure has been imple¬ 
mented. The chosen preconditioner is based on the ISF 
Poisson solver for the standard Poisson equation, which 
can handle all common boundary conditions exactly. The 
code requires a small amount of memory and is very fast 
and in addition also highly parallelized. We have shown 
that coupled with BigDFT, our method can correctly re¬ 
produce electrostatic solvation energies of a set of small 
neutral organic molecules. Effects of grid resolution and 
boundary conditions on the electrostatic solvation ener¬ 
gies have been reported. The whole library, will be re¬ 
leased as an independent program suitable for integration 
in other electronic structure codes. A GPU-accelerated 
version of this software package is also in preparation, 
following the guidelines indicated in Ref. [ 37]. 
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